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ABSTRACT 



We study the properties of galaxy super- winds driven in major galaxy mergers, using high-resolution 
(~ 1 pc) hydrodynamic simulations with detailed models for stellar feedback that can follow the formation 
and destruction of individual GMCs and star clusters and self-consistently follow the generation of winds. 
The models include molecular cooling, star formation only at extremely high densities inside GMCs, and 
gas recycling and stellar feedback from SNe (Type I & II), stellar winds, and radiation pressure (both 
in UV/optical and re-processed IR). We study mergers of a range of systems from SMC-like dwarfs and 
Milky Way analogues to z ~ 2 starburst disks. In all cases, multi-phase super- winds are generated in all 
passages, with peak outflow rates from ~ 10 — 1000 Mq yr _1 . However, the wind mass-loading efficiency 
(wind mass loss rate divided by the star formation rate) is similar to that in the isolated galaxy counterparts 
of each merger: it depends more on global galaxy properties (mass, size, and escape velocity) than on the 
dynamical state of the galaxy or parameters of the merger. The winds tend to be bipolar or even unipolar, 
but multiple overlapping "events" generated by passages lead to complex morphologies with overlapping 
and differently-oriented bubbles/shells at a range of radii. The winds have complex velocity and phase 
structure, with material at a wide range of outflow speeds up to ~ 1000km s _1 , and a mix of molecular, 
ionized, and hot gas that depends on the galaxy properties; these different phases are connected to different 
feedback mechanisms. These simulations resolve a long-standing problem in "sub-grid" models: models 
with simple wind prescriptions typically find that they dramatically suppress merger-induced starbursts, 
often making it impossible to form ULIRGs even in major mergers and eliminating any significant merger- 
induced SF enhancement in gas-rich galaxies. But despite large mass-loading factors up to ~ 10 — 20 in 
the winds, the peak SFRs in the simulated starbursts here are as high or higher than those in "no wind" 
simulations. The acceleration of winds does not act equally on all material, so cold, dense gas which is 
strongly torqued can still lose angular momentum and contribute to a nuclear starburst while blowing out 
diffuse, low-density gas that would not have participated in the starburst in the first place. A significant 
fraction of the wind material is not entirely unbound, and falls back into the disk at later times after the 
merger, leading to higher post- starburst SFRs in the presence of stellar feedback (compared to models 
with no winds at all). This may in turn require strong AGN feedback to explain galaxy "quenching." 

Key words: galaxies: formation — galaxies: evolution — galaxies: active — star formation: general — 
cosmology: theory 



1 INTRODUCTION 

Feedback from massive stars is critical to the evolution of galaxies. 
In cosmological models of galaxy evolution lacking strong stellar 
feedback, gas cools and turns into stars rapidly, producing galaxies 
with order-of-magnitude larger stellar masses than those observed 



(e.g. Katz et al.|1996||Somerville & Primack 1999 ; 


Coleetal.|2000; 


Springel & Hernquist||2003b| |Keres et al.|2009a 


and references 



therein). Simply 'forcing' inefficient star formation does not re- 
solve this; the amount of baryons in real galactic disks is much 
lower than the amount of cool gas in disks found in cosmologi- 



cal simulations, especially in low-mass galaxies (White & Frenk 
|1991[ for a recent review see Keres et al. 2009b). From the ob- 
served mass-metallicity relation and enrichment of the IGM, it is 
also clear that the baryons must have at one point entered galaxy 
halos and disks and been enriched, then were ejected (Tremonti 
[eTaT1l2004l |Erb et al ||2006| |Aguirre et al.| [200T] |Pettini et"alj 
2003, Songaila 2005, Martin et al.||2010| )7Galactic super- winds 
are therefore implied, with large mass-loading factors of several 
times the SFR that are required in cosmological simulations to re- 
produce these observations (e.g. Oppenheimer & Dave 2006). Such 
mass-loading factors are also observationally inferred in many lo- 
cal galaxies and massive star-forming regions at z ~ 2 — 3 (Martin 
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T999l|2006l[Heckman et al.|2000l|Newman et al.|20T2l|Sato et all 



2009l|Chen et~aT 2010, Steidel et al. 2010, Coil et al.|20lT) . 



Until recently, however, numerical simulations have generally 
been unable to produce, from an a priori model, winds with large 
mass-loading factors (as well as a plausible scaling of wind mass- 
loading with galaxy mass or other properties); this is especially true 
of models which include only thermal or "kinetic" feedback via 
supernovae, which is very inefficient in the dense regions where 



star formation occurs (see e.g. 


Guo et al.|2010||Powell et al.|2011| 


|Brook et al.||2011 ; Nagamine 


2010 and references therein). This 



should not be surprising: feedback processes other than supernovae 
are critical for suppressing star formation in dense gas; these in- 
clude protostellar jets, HII photoionization, stellar winds, and radi- 
ation pressure from young stars. Including these mechanisms self- 
consistently maintains a reasonable fraction of the ISM at densi- 
ties where the thermal heating from supernovae has a larger effect; 
moreover there are many regimes where these mechanisms can di- 
rectly drive winds, independent of and with greater mass loading 
than supernovae. 

This conclusion implies that (not surprisingly) an accurate 
treatment of galactic winds requires a more realistic treatment of 
the stellar feedback processes that maintain the multi-phase struc- 
ture of the ISM of galaxies. Motivated by these problems, in |Hop-| 
|kins et al.| ( [20TTbt (Paper I) and |Hopkins"eTaT] pOT2d] ) (Paper II), 
we developed a new set of numerical models to follow feedback on 
small scales in GMCs and star-forming regions, in simulations with 
pc-scale resolution^These simulations include the momentum im- 
parted locally (on sub-GMC scales) from stellar radiation pres- 
sure, radiation pressure on larger scales via the light that escapes 
star-forming regions, HII photoionization heating, as well as the 
heating, momentum deposition, and mass loss by SNe (Type-I and 
Type-II) and stellar winds (O star and AGB). The feedback is tied 
to the young stars, with the energetics and time-dependence taken 
directly from stellar evolution models. Our models also include 
cooling to temperatures < 100K, and a treatment of the molecu- 
lar/atomic transition in gas and its effect on star formation (follow- 
ing [Krumholz & Gnedin 2011). We showed that these feedback 
mechanisms produce a quasi- steady ISM in which giant molecular 
clouds form and disperse rapidly, after turning just a few percent 
of their mass into stars. This leads to an ISM with phase structure, 
turbulent velocity dispersions, scale heights, and GMC properties 
(mass functions, sizes, scaling laws) in reasonable agreement with 
observations. 

In Hopkins ~et al.| ( [2012c] ) (Paper III), we show that these same 
models of stellar feedback predict the elusive winds invoked in 
almost all galaxy formation models; the combination of multiple 
feedback mechanisms is critical to give rise to massive, multi-phase 
winds having a broad distribution of velocities, with material both 
stirred in local fountains and unbound from the disk. 

However, in Paper III we examine only idealized isolated disk 
galaxies. Although this is probably representative of much of a 
galaxy's lifetime, a great deal of observational study has focused 
on winds in "starburst" galaxies, often in interacting or merging 
systems. Indeed, a wide range of phenomena indicate that gas- 
rich mergers are important to galaxy formation and star forma- 
tion. These systems dominate the most intense starburst popula- 
tions: ULIRGs at low redshift ([Joseph & Wright|1985||Sanders &[ 



Mirabel 1996, Evans etal.|2009 



and Hyper-LIRGs and bright sub- 



millimeter galaxies at high redshifts Pap ovich et al. 12005] [You nger 



|et al.|2008b||Tacconi et al.|2006[|Sctunnerer et al.|2008||Chap man 



et al. 2009, Tacconi et al. 2008 ). They are powered by compact 
concentrations of gas at enormously high densities (Scoville et al. 
[T9861 |Sargent et aT][T987] ), which provides material to fuel BH 
growth and boost the concentration and central phase-space den- 
sity of merging disks to match those of ellipticals (Hernquist et al. 
[T9^|Robertson et al.|20061|Hopkins et al.|2008b| ). Various stucl 
ies have shown that the mass involved in these starburst events is 
critical for explaining the relations between spirals, mergers, and 
ellipticals, and has a dramatic impact on the properties of merger 
remnants (e.g., Lake & Dressle r|1986[|Doyon et al.|19"94{ [Sh ier & 
Fischer|1998HJames et al.| 19991 |Genzel et al.|200l||Tacconi et al.| 



2002, |Dasyra"eTaL| 2006 2007, Rothbe rg & Josep h 2004 2006, 



1 Movies of these simulations are available at jhttp : / /www .tapir . | 
[caltech . edu/ ~phopkins/Site/Movies_sbw_mgr . html 



Hopk ins et al.|2009a|b| >. 

With central densities as large as ~ 1000 times those in Milky 
Way giant molecular clouds (GMCs), these systems also provide a 
laboratory for studying star formation, the ISM, and the generation 
of galactic winds under the most extreme conditions. In Hopkins 
|et al.| ( [2012a| > (Paper IV), we therefore extend the models from Pa- 
per I-Paper III to include major galaxy mergers. We showed there 
that the same feedback mechanisms can explain the self-regulation 
of starbursts and extension of the Kennicutt-Schmidt relation to 
the highest gas surface densities observed. We also show how this 
controls the star formation rates and their spatial distributions, the 
formation of clusters, and the formation and destruction of giant 
molecular clouds in the ISM. In this paper, we further investigate 
the phase structure and generation of galactic superwinds in these 
models, and how they relate to merger dynamics and star formation 
histories. 

2 METHODS 

The simulation techniques and galaxy models used here are de- 
scribed in detail in Paper I (Sec. 2 & Tables 1-3) and Paper II 
(Sec. 2). The details of the merger models are given in Paper IV. We 
briefly summarize the most important properties here. The simula- 
tions were performed with the (heavily modified) parallel TreeSPH 
code GADGET-3 (Springel 2005 ), in its fully conservative formula- 
tion (Springel & Hernquist 2002 ). They include stars, dark matter, 
and gas, with cooling, shocks, star formation, and stellar feedback. 

2.1 Disk Models 

We implement the ISM prescription in four distinct initial disk 
models spanning a range of galaxy types. Each has a bulge, stel- 
lar and gaseous disk, halo, and central BH (although to isolate the 
role of stellar feedback, models for BH growth and feedback are 
disabled). At our standard resolution, each model has ^0.3 — 1 x 
10 8 total particles, giving particle masses of 500 — 1000 Mq and 
1 — 5 pc smoothing lengths. A couple of ultra-high resolution runs 
(of isolated versions of these disks) for convergence tests employ 
w 10 9 particles with sub-pc resolution|^The disk models include: 
(1) SMC: an SMC-like dwarf, with baryonic mass Mb ar = 



2 These tests are described in Paper I and Paper II, and used to check con- 
vergence in small-scale ISM properties. We have also run every simulation 
described in this paper with 10 times fewer particles (2 times larger soft- 
ening); although the small-scale properties such as the structure of indi- 
vidual star clusters differ, all integral properties plotted in this paper (more 
global values in general) are similar. Additional numerical tests of the SPH 
method, relevant primarily for the wind phase structure in the extended halo, 
are presented in Appendix [A| 
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Figure 1. Morphology of the gas in a standard simulation (high-resolution 
with all explicit feedback mechanisms included) of a merger of the HiZ disk 
model: a massive, z ~ 2 — 4 starburst disk merger. The time is near apocen- 
ter after first passage. This color projection emphasizes the ionized and hot 
gas. Brightness encodes projected gas density (increasing with density; log- 
arithmically scaled with a w 4dex stretch); color encodes gas temperature 
with the blue/white material being T < 1000K cold atomic/molecular gasP 
pink - 10 4 - 10 5 K warm ionized gas, and yellow > 10 6 K hot gas. Violent 
outflows are present, emerging both from the complexes and the system as 
a whole, driven by the massive starburst. The volume-filling "hot" compo- 
nent is now visible, with multiple "bubbles" driven by separate local events. 



t = 0.6 




Figure 2. The stars at the same time as Fig. [T] The image is a mock ugr 
(SDSS-band) composite, with the spectrum of all stars calculated from 
their known age and metallicity, and dust extinction/reddening accounted 
for from the line-of-sight dust mass. The brightness follows a logarithmic 
scale with a stretch of « 2 dex. Young star clusters are visible throughout the 
system as bright white pixels. The nuclei contain most of the star formation, 
evident in their saturated brightness. Fine structure in the dust gives rise to 
complicated filaments, dust lanes, and patchy obscuration of star-forming 
regions. A few super-star clusters are apparent as the brightest young stellar 
concentrations. 



8.9 x lO ^M^ and halo ma ss Mhaio = 2 x 10 10 Mq (concentration 
c — 15), a Hernquist ( 1990) profile bulge with a mass mb = 10 7 Mq, 
and exponential stellar {m d = 1.3 x 10 8 Mq) and gas disks (m g — 
7.5 x 10 8 Mq) with scale-lengths h d — 0.7 and h g = 2.1kpc, re- 
spectively. The initial stellar scale-height is zo = 140 pc and the stel- 
lar disk is initialized such that the Toomre Q—\ everywhere. The 
gas and stars are initialized with uniform metallicity Z = O.IZq. 

(2) MW: a MW-like galaxy, with halo and bary- 
onic properties of (Mhaio, c) — (1.6 x 10 12 Mq, 12) and 
(M har ,m b ,m d ,m g ) = (7.1,1.5,4.7,0.9) x 1O 1O M , Z = Z , 
and scale-lengths (h d ,h g ,zo) = (3.0, 6.0, 0.3) kpc (note the gas 
disk is more extended than the stellar disk, giving a gas fraction 
« 10% inside the solar circle). 

(3) Sbc: a LIRG-like galaxy (i.e. a more gas-rich spiral than is 
characteristic of those observed at low redshifts) with (Mhaio, c) — 
(1.5 X 10 n MQ, 11), (Mbar, m&, m d , m g ) — (10.5, 1.0, 4.0, 5.5) x 
10 9 M , Z = O.3Z , and (h d , h 8 , zo) = (1.3, 2.6, 0.13) kpc. 

(4) HiZ: a high-redshift massive starburst disk, chosen to 
match the properties of the observed non-merging but rapidly star- 
forming SMG population, with (M ha io, c) = (1.4 x 10 12 Mq, 3.5) 
and a virial radius appropriately rescaled for a halo at z = 2 rather 
than z = 0, (M ba r, m b , m d , m g ) (10.7, 0.7, 3, 7) x 1O 1O M , Z 
0.5Zq, and (h d , h g , zo) = (1.6, 3.2, 0.32) kpc. 

2.2 Merger Parameters 

The purpose of this paper is to study the generation of winds in 
mergers by stellar feedback, not to present a systematic study of all 
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merger parameters - extensive studies of this nature can be found in 
the lit erature ([Barnes & Hernquist|1992l|Barnes|2002[|Naab et al. 
|T999l|CoTet al.|2006a||2008||Robertson et al.|2006|| Younger et al. 



2008a, di Matteo et al. 2007, Burkert et al. 2008, Hoffman et al. 



2009 , Jess eit et al.|2 009 ). We therefore focus on a small but repre- 
sentative subset of possible mergers. 

We consider equal-mass mergers (merging identical copies of 
galaxies (l)-(4)). The pairs are placed on parabolic orbits (mo- 
tivated by cosmological simulations; see Benson 2005] |Khoch-| 
|far & Bur kert 2006) with the spin axis of each disk specified 
by and 4> in spherical coordinates. We choose two represen- 
tative orientations from the which have been studied in a num- 
ber of othe r works (|Cox et al.||2006a[ |Barnes|p"988[ IHern quist 
^992l|T993l|Naab & Burkert|2003| ). The first (orbit e in lCox et al. 
2006a] is near-prograde with (0 1? <j> u 2 , fa) = (30, 60, -30, 45); 
the second (orbit f) is near-retrograde (or polar-retrograde) with 
(0i, 0i, 02, fa) = (60, 60, 150, 0). This spans both a strong res- 
onant interaction (prograde) and weak, out-of-resonance (retro- 
grade) interaction (e.g. jD'Onghia et al. 2009 2010). We choose 
these rather than perfectly prograde (0i, fa, 02, fa) — (0, 0, 0, 0) 
and retrograde (0i, fa, 62, fa) = (180, 0, 180, 0) orbits because 
there are known pathological behaviors associated with the narrow 
phase space of perfectly resonant orbits. 

We will show that, for the relevant properties of stellar winds, 
there is little difference between these orbits which bracket the 
range from most to least violent encounters; we therefore expect 
the properties examined here to be robust over a wide range of con- 
figurations. 
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HiZ 




Figure 3. Images of the gas (color-coded as Fig.^ in the starburst-driven superwinds in a our simulations. We show just one example of each merger, but 
the e and f orbits are similar in each case. For each, we show the image at a fixed time near the starburst and at various spatial scales. A massive wind is 
plainly evident; the winds are multi-phase with volume-filling hot gas and ejected streams/shells of cold gas. We caution that the formation of the small, 
marginally-resolved isolated "blobs" of cold gas within the hot background outflow (not the warm/cold shells or GMCs within the disk) is quite sensitive to 
the numerical details of the simulations!^ A large fraction is unbound and escapes the galaxy halo. In the MW case, the gas-poor nature of the merger means 
the wind is almost entirely "hot"; in the Sbc and HiZ cases a much larger fraction of ejected material is warm/cool; in the SMC case there is a broad mix of 
warm gas and hot gas ejected. The shells and features in the diffuse gas arise from multiple bursty episodes shocking as they "catch up" to one another. 



2.3 Star Formation and Feedback 

The most important physics in these simulations is the model of 
stellar feedback. We include feedback from a variety of mecha- 
nisms, each of which we briefly describe below. More details about 
our implementations of this physics are given in Paper I and Paper 
II. We use a|Kroupa (2002]) initial mass function (IMF) through- 
out and use STARBURST99 ( |Leitherer etal.|1999| > to calculate the 
stellar luminosity, mass return from stellar winds, supernova rate, 
etc. as a function of the age and metallicity of each star particle. 



(1) Local Momentum Deposition from Radiation Pressure, 
Supernovae, & Stellar Winds: In Paper I, we present the radiation 
pressure aspect of this model for feedback from young star clusters 
in detail. At each timestep, gas particles identify the nearest density 
peak representing the center of the nearest star-forming "clump" or 
GMC-analog. We calculate the total luminosity of the star particles 
inside the sphere defined by the distance from the center of this 
star-forming region to the gas particle of interest; the incident flux 
on the gas is then determined assuming that the local star forming 
region is optically thick to the UV radiation. 
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The rate of momentum deposition from radiation pressure is 
then P r ad ~ (1 + tir) Li nc ident/c where the term 1 + tir accounts for 
the fact that most of the initial optical/UV radiation is absorbed and 
re-radiated in the IR; tir = S ga s ftiR is the optical depth in the IR, 
which allows for the fact that the momentum is boosted by multiple 
scatterings in optically thick regions. Here Hg as is calculated self- 
consistently as the average surface density of the identified clump, 
with kir w 5 (Z/Zq) g _1 cm 2 approximately constant over the rel- 
evant physical range of dust temperatures. The imparted accelera- 
tion is directed along the flux vector. In Paper I we discuss numer- 
ous technical aspects of this implementation - such as the effects of 
resolution, photon leakage, and how the momentum is discretized 
- and show that essentially all our conclusions are robust to uncer- 
tainties in these choices. 

The direct momentum of SNe ejecta and stellar winds PsNe 
and P w are similarly tabulated from STARBURST99 and injected 
as an appropriate function of age and metallicity to the gas within 
a smoothing length of each star. This source of turbulent energy is 
almost always smaller than that due to radiation pressure discussed 
above. In some cases, however, in particular in dwarf galaxies, the 
work done by bubbles of gas shock-heated by supernovae and/or 
stellar winds is dynamically important; this is discussed below. 

(2) Supernova and Stellar Wind Shock-Heating: The gas 
shocked by supernovae and stellar winds can be heated to high tem- 
peratures, generating bubbles and filaments of hot gas. We tabulate 
the Type-I and Type-II SNe rates from Man nucci et al.| ( |2006] ) and 
STARBURST99, respectively, as a function of age and metallicity 
for all star particles, and stochastically determine at each timestep 
if a SN occurs. The SNe are resolved discretely in time (as opposed 
to continuous energy injection). For each SN, the appropriate ther- 
mal energy is injected into the gas within a smoothing length of 
the star particle. Similarly, stellar winds are assumed to shock lo- 
cally and so we inject the appropriate tabulated mechanical power 
L(t, Z) as a continuous function of age and metallicity into the gas 
within a smoothing length of the star particles. The specific energy 
of these stellar winds is large for young stellar populations in which 
fast winds from massive stars dominate, but declines rapidly at later 
times when slower, AGB winds dominate. 

(3) Gas Recycling: Gas mass is returned continuously to the 
ISM from stellar evolution, at a rate tabulated from SNe and stellar 
mass loss in STARBURST99. The integrated mass fraction recy- 
cled is ~ 0.3. 

(4) Photo-Heating of HII Regions: For each star particle, we 
tabulate the rate of production of ionizing photons; starting from the 
nearest gas particle and moving radially outwards, we then ionize 
each non-ionized gas particle (using the gas and stellar properties 
to determine the photon production rate needed to maintain the par- 
ticle as fully ionized). Gas which is ionized is immediately heated 
to ~ 10 4 K, unless is already above this temperature; moreover, the 
gas is not allowed to cool below 10 4 K until it is no longer in an 
HII region. This method allows for overlapping, non-spherical HII 
regions that can extend to radii ^kpc. 

(5) Long-Range Radiation Pressure: Radiation pressure 
from photons absorbed in the immediate vicinity of stars is cap- 
tured in mechanism (1). However, photons that escape these re- 
gions can still be absorbed at larger radii. For each star particle, we 
construct the intrinsic SED (L u ) as a function of age and metallic- 
ity; we then use the local density and density gradients to estimate 
the integrated column density and attenuation of the SED using 

p(h sm \ + | Vlnp| l ), where /z sm i is the smoothing 
length and k v is the frequency-dependent opacity (assuming dust 
opacities that scale with metallicity, as in (1) above). The resulting 



"escaped" SED gives a frequency-dependent flux F u that is prop- 
agated to large distances. We construct a force tree for this long 
range force in an identical fashion to the gravity tree, since after 
attenuating the flux near the star particle, the stellar flux is assumed 
to decrease oc r -2 . Each gas particle then sees an incident net flux 
vector F z ^, integrated over all stars in the galaxy. 

Extensive numerical tests of the feedback models are pre- 
sented in Paper I, Paper II, and Paper III. 

In these models, gas follows an atomic cooling curve with ad- 
ditional fine- structure cooling to < 100K, with no cooling floor 
imposed. 

Star formation is allowed only in dense, mo lecular, self- 
gravitating regions above n > 1000 cm -3 . We follow Krumholz & 
Gnedin (201 1) to calculate the molecular fraction fu 2 in dense gas 
as a function of local column density and metallicity, and allow SF 
only from molecular gas. Following Hopkins et al. ( 201 la), we also 
restrict star formation to gas which is locally self-gravitating, i.e. 
has a = 6v 2 6r/Gm gas (< Sr) (| V • v| 2 + | V x \\ 2 )/(2Gp) < 1 
on the smallest available scale (Sr being our force softening or 
smoothing length). This forms stars at a rate p* = p mo \/tn (i-e. 
100% efficiency per free-fall time); however the average efficiency 
is much lower (~ 1 %) because of feedback. In Paper I and Paper 
II we show the galaxy structure and SFR are basically independent 
of the small-scale SF law, density threshold (provided it is high), 
and treatment of molecular chemistry, because they are feedback- 
regulated. 

In Paper I and Paper II we show that the galaxy structure and 
SFR are basically independent of the small-scale SF law, density 
threshold (provided it is high), and treatment of molecular chem- 
istry, because they are feedback-regulated. We have, for exam- 
ple, re-run several of the models here with a more traditional pre- 
scription where star formation is restricted to gas above a density 
n > 1000cm -3 , with a rate = 0.015 p/t$ corresponding to ob- 
served average efficiencies "enforced." This has almost no effect 
on global properties or the galaxy- average SFRs, but does tend to 
"smear out" the SFR in dense regions such as the merging nuclei, 
where all the gas is well above this threshold. 

3 OUTFLOW MORPHOLOGIES 



Figures |1|2| show the gas and stars in a representative stage of the 
merger simulations, when all feedback mechanisms are present. 
The gas maps show the projected gas density (intensity) and tem- 
perature (color, with blue representing cold atomic/molecular gas 
at T < 1000K, pink representing the warm ionized gas at ~ 
10 4 — 10 5 K, and yellow representing the hot, X-ray emitting gas 
at > lO 6 K)0The stellar maps show a mock three-color observed 
image, specifically au/g/r composite]^] 

Outflows are evident; we will quantify their properties and 
phase distribution in detail below. Briefly, the simulated ISM is a 
supersonically turbulent medium (see Hopkins 2012b c) in which 

3 The projected temperatures are logarithmically-averaged and surface- 
density weighted, so reflect the temperature of most of the line-of-sight gas 
mass, rather than the temperature that contains most of the thermal energy. 

4 The stellar luminosity in each band is calculated from each star particle 
according to the STARBURST99 model given its age, mass, and metallic- 
ity (and smoothed over the appropriate kernel). We then attenuate the stars 
following the method of Hopkins et al. (2005): we calculate the total dust 
column (from the simulated gas) along the line-of-sight to each star parti- 
cle for the chosen viewing angle (assuming a constant dust-to -metals ratio, 
i.e. dust-to-gas equal to the MW value times Z/Zq), and apply a MW-like 
extinction and reddening curve (as tabulated in Pei 1992). 



© 0000 RAS, MNRAS 000, 000-000 




Figure 4. Galactic wind thermal+metal line emission morphology, as a proxy for X-ray emission (though we caution this includes gas with a broad range 
of temperatures). The maps show the images as Fig. [3] but in a single-color scale where intensity here encodes the projected bremsstrahlung emissivity plus 
metal cooling luminosity (again we caution that the small clumps at large radii may be artificiaPj). Since this weights the volume-filling hot gas, the different 
bubbles and shells are more clear. At each scale, clear changes in the orientation of the older/larger bubbles are also evident. 



cold giant molecular clouds (GMCs) continuously form and are 
dispersed by feedback after turning a few percent of their mass 
into stars (see Paper II, Paper IV). There is volume-filling hot gas 
(heated by SNe and O-star winds), which vents an outflow compo- 
nent; occasionally the early stages are evident as "bubbles" break- 
ing out of the warm phase gas. Warm/cold gas is mixed through- 
out the outflowlj especially near the galaxy (before they have had 



5 We caution that the very small, cold "clumps" at large radii in the wind 



time to mix with the more diffuse material); these are sometimes 



(distinct from the resolved molecular clouds in the disk and the coherent 
warm/cold shells and filaments in outflow) are only marginally resolved 
and are probably an artifact of the numerical method used, in concert with 
the fact that there is no cosmic ionizing background present in these simu- 
lations to heat gas far from the galaxy. In Appendix |A| we re-run a limited 
subset of our simulations with an alternative formulation of SPH designed 
to treat contact discontinuities more accurately, and include UV background 
heating, and show that these largely disappear. We will focus in this paper 
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entrained by the hot gas but more often directly accelerated by ra- 
diation pressure. Qualitatively, this is true in isolated galaxies as 
well (Paper III). 

Figure [3] shows images of the gas in these winds at progres- 
sively larger radii. In the images, successive generations of "bursts" 
and strong outflows are evident as overlapping shells (many caused 
by shocks as different ejection events "catch up"); these are each 
broadly associated with a galaxy passage. The diffuse hot gas, 
being volume-filling, has a near-unity covering factor, but it is 
still clearly organized into various "bubbles." Each of these has 
a broadly bimodal morphology set by much of the hot gas blow- 
ing out perpendicular to the disk along the path of least resistance. 
Since the disk orientations change during the merger owing to grav- 
itational torques, the successive bursts have different orientations. 
At the largest radii, the diffuse gas cools adiabatically; this would 
not necessarily occur with a realistic IGM present into which the 
hot gas could propagate. The warm and cool gas also has a signif- 
icant (albeit smaller) covering factor even at > lOOkpc (except in 
the MW model). However, the morphology is much less smooth - 
this gas is primarily in filaments, and shells. Although the crude 
average distribution of these is similar to the hot gas, there is much 
larger line-of- sight variation. Some (but not all) of this is material 
is accelerated by radiation pressure, and so tends to reach some- 
what lower velocities than the hottest pressure-accelerated gas, and 
therefore the density falls off proportionally more rapidly at very 
large radii (although, being more dense, this material may be able to 
propagate ballistically into the IGM where lower-density, volume- 
filling material would be halted by the ambient pressure). In any 
case, as we discuss in detail in Paper III and §[8] the outflows should 
not be taken too literally at the largest radii: we do not include ini- 
tial gaseous halos or a realistic IGM into which the wind should 
propagate, so the winds here expand unimpeded beyond the halo. 

X-ray observations provide a strong probe of the hot phase 
of the galactic winds; Figure [4] therefore shows the same images, 
but now in their approximate X-ray properties. For convenience, 
rather than make a detailed mock observation corresponding to 
a given instrument, sensitivity, redshift, and energy range, we in- 
stead quantify the approximate X-ray emission with the sum of 
the thermal bremsstrahlung emission (emissivity per unit volume 
ux oc TgJ s 2 n e rii, in terms of the gas temperature r gas , and elec- 
tron/ion number densities n e /rii), and the metal-cooling luminos- 
ity (using the compiled tables in Wiers ma et al.| |2009 as a func- 
tion of Hi, n e , Z, and T gSLS , assuming solar abundance ratios and the 
z — ionizing background). Note that these can include signifi- 
cant contributions from low-temperature gas so this need not re- 
fer specifically to X-ray observations. Broadly, the morphology is 
similar, but with the hot, low-density material highlighted, the var- 
ious bubbles and shells are more obvious. It is also clear, as we go 
to larger scales, that the larger/older bubbles have distinct orien- 
tations, corresponding to the orientation of the galaxies at earlier 
merger stages. 

4 PHASE STRUCTURE: HOT X-RAY HALOS AND COLD 
MOLECULAR GAS 

In Paper II, we discuss in great detail the phase structure and den- 
sity distribution of the ISM; here, we examine whether the same 
results obtain in mergers. Fig.[5]plots the density PDF of the ISM, 
i.e. mass per logarithmic interval in n. For clarity, we just show 

only on results robust to our resolution and numerical method, but direct the 
interested reader to AppendixlAland Hopkins (2012a| for more details. 



the f models, but the e models are extremely similar. This covers 
a very wide dynamic range and is not directly comparable to the 
typical "ISM density distribution" from galaxy studies. To see this, 
we divide the gas into three categories: the "star-forming disk(s)" 
(gas within R^, within one exponential "scale height" defined with 
respect to the angular momentum plane, and with outflow veloc- 
ity v r < 100km s _1 ); the wind/outflow (defined below as un-bound 
gas with large outflow velocity), and extended disk+halo gas (the 
remaining gas; recall there is no initial extended gaseous halo in 
these simulations). Unsurprisingly, the SF "disks" include most of 
the dense gas, the winds include the least dense material (much of 
it out at or past the virial radius), and the "halo" is intermediate. We 
note again that since there is no IGM, escaped material can reach 
arbitrarily low densities. We also compare the distribution of den- 
sity (weighted by mass or volume) in different temperature phases 
corresponding to cold/warm ionized/hot X-ray emitting gas within 
the star forming regions/disks (the exact temperature cuts are ar- 
bitrary but the qualitative comparison does not change if we shift 
them by moderate amounts). Here, we recover a qualitatively sim- 
ilar result to the isolated disks: the high-density gas (> 1 cm -3 ) is 
predominantly cold and contains most or ^half the mass but has 
a small volume-filling factor of a couple percent; the intermediate- 
density gas (0.01 — 1 cm -3 ) is primarily warm and has both a signif- 
icant fraction of the mass (~ 30 — 50%) and sizeable filling factors; 
the low-density gas (< 0.01 cm -3 ) is primarily hot and has order- 
unity volume filling factors. Each component is crudely (but not ex- 
actly) log-normal. As discussed in Paper II, the turbulent pressure 
inside GMCs is much larger than the background pressure (they are 
marginally self-gravitating, rather than pressure-confined). 

Fig. [6] shows how this compares (over the course of the 
merger) to the isolated disks. We show one example but the re- 
sults are similar in each case. We plot the density distribution by 
temperature phase for all gas, averaged over the run of the isolated 
disk simulation (in which it reaches a steady- state, so is nearly 
time-independent), averaged over the merger duration, and at the 
snapshot with the peak SFR (near final coalescence) in the merger. 
Averaged over the merger, there is little difference (some material 
is at lower densities simply because the winds have more time to 
"escape," and slightly more material is at high densities). This re- 
flects the well-known fact that in merger simulations, most of the 
time (and most of the star formation) is contributed by the "iso- 
lated mode," namely the separate SF in the two disks as they orbit 
between passages, rather than the merger-induced "burst" on top 
of this (which only dominates the central ^kpc seen in Paper IV; 
see also |Cox et al.|2008| [Hopkins & Hernquist|20T0] >. During the 
peak starburst, the results are still qualitatively similar (especially 
for the warm/hot gas); the peak density associated with each phase 
is also the same. The main difference is that a fraction of the gas 
funneled into the galaxy centers is pushed to much higher densi- 
ties n ^> 10 4 cm -3 where it rapidly turns into stars. Although we 
stress that this is still not most of the gas (even the dense molecu- 
lar gas). This is very different from models with weak/no feedback, 
which see catastrophic runaway to arbitrarily high densities n ~ 10 6 
for most of the star-forming gas, in contrast to observations which 
show that most of the gas in GMCs and other dense regions is at 
modest, non star-forming densities (e.g. Williams & McKee 1997, 
|Evans|1999l and references therein). 

The increase in gas at the highest densities during the star- 
burst would be evident in dense molecular traces, as discussed 
in |Narayanan et al.| ( |2006| ); |Hopkins et al.| ( |2012b| ). Specifically, 
if we adopt the simple conversions therein from mass above ~ 
10 4 cm -3 to HCN luminosity and mass above ~ 100 cm -3 to CO(l- 
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Figure 5. Density distribution of different ISM phases, for different galaxy disk models (Sbc, MW, and SMC), averaged over the duration of the merger. Top: 
Mass-weighted density PDF (dm gas /dlog«), i.e. the mass fraction per logarithmic interval in density n. We show the distribution for all gas in the simulation 
(black), the gas approximately within the (multi-phase) star-forming disks (blue), the gas in the extended, ionized disks and halo (green), and wind/outflow 
material (orange). These trace the material at low, intermediate, and high densities, respectively (as expected). Each density PDF has a very broad density 
distribution. Middle: Mass-weighted density PDF within the "star-forming disks". We show the density PDF for all of the gas in the region (black), the 
"cold" phase (purple; T < 2000 K), the "warm ionized" phase (cyan; 2000 < T < 4 x 10 5 K), and the "hot diffuse" phase (red; T > 4 x 10 5 K). Most of the 
mass is in the cold phase (which dominates at high densities), but with a comparable contribution from the warm medium (the hot phase contributing a few 
percent). The multi-modal nature of the total density PDF is a consequence of the strong phase separation. Note that since these simulations do not include the 
cosmic ionizing backgro und, we truncate the cold-phase distribution at densities below which post-processing calculations suggest it would be photo-ionized 
(~ 0.01cm -3 ; see Faucher-Giguer e et al. 2008); such gas contributes negligibly to the total at very low densities, in any case. Bottom: Volume-weighted 
density PDF (dVg as /dlogft) for the star forming disk. The hot diffuse phase dominates, with a moderate volume filling fraction for the warm phase, and a 
small (~ 1—5%) volume filling fraction of cold molecular clouds in the disk. 



0) luminosity, we estimate that the ratio Lhcn/^co(i-o) should 
increase from « 0.02 in the isolated case (presented therein) to 
~ 0.15 — 0.30 in the peak of the merger. This is consistent with 
what is seen in real local ULIRGs, essentially all of which are late- 
stage major mergers, in the compilatio n of|Gao & Solomon| ( |2004| ); 
|Narayanan etaTl ( [2008l ); |Juneau et al.|f2009| >. 

In Fig^] we show the phase distribution of the winds, sim- 
ilar to Fig. [5] These are discussed in detail for the isolated cases 
in Paper III. We compare the distribution of temperatures weighted 
by their contribution to thermal bremsstrahlung plus metal cooling 
emission (as described above) for all the gas, in different outflow 
velocity intervals. We also plot the velocity distribution of all gas 
(mass per radial outflow velocity v r ); the "peak" near v r — be- 
ing the non-wind material (with the wind evident in the large tail). 



And we show the density distribution of the wind material specifi- 
cally, in the style of Fig. [5] The velocity distributions are wide (dis- 
cussed below). The inhomogeneous morphologies of the winds are 
reflected in the broad phase and density distribution (though recall 
the numerical caveats regarding the clumpy structure of the out- 
flows). The extremely low-density wind material is a consequence 
of our not including a full IGM into which the winds can prop- 
agate; this also causes the winds to cool adiabatically (hence the 
secondary "bump" of warm/cold gas at extremely low densities). 
However it is fairly generic that the cold/warm/hot material domi- 
nates the wind at high/intermediate/low densities, respectively. It is 
especially worth noting that the winds can include some cold gas 
at densities ~ 1 — 100 cm -3 (this is primarily material still near or 
within the disk, being accelerated outwards; as it escapes, the gas 
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Figure 6. Density distribution of different ISM phases (as Fig. [5j, for all 
the gas in the simulation in one example (the Sbc e case, but the results are 
qualitatively similar for each model). Top: The average, steady-state distri- 
bution for the isolated disk (see Paper II). Middle: Distribution averaged 
over the duration of the merger. Bottom: Distribution at the snapshot where 
the merger-induced SFR is maximized. The average phase distribution over 
the entire merger is similar to (just slightly more broad than) the isolated 
disk counterpart, reflecting the fact that most of the SF is in the separate 
"quiescent mode" in the disks rather than the merger-induced burst. Unlike 
a no-feedback model, there is not a runaway pileup of gas at the highest 
densities. During the peak of activity, more gas is channeled to high densi- 
ties > 10 4 cm 5 , but otherwise the distributions (and characteristic density 
of each temperature phase) are similar. This would be observable in e.g. the 
ratio of £hcn/£cq(i-0) which would increase from « 0.02 in the isolated 
case (see Hopkins et al. 2012b) to « 0.2 — 0.3 in the peak of the merger. 



expands and can easily be heated) P In the gas-rich mergers, the 
contribution to the wind mass from cold/warm/hot phases are com- 
parable, although in the dwarf systems (SMC and Sbc) the very 
low-density "warm" component was mostly "hot" when originally 
ejected; in the gas-poor MW-like case, there is relatively little ma- 
terial to "entrain" and torques are able to efficiently force most of 
the dense gas into a starburst, so the outflow is much more strongly 
dominated by hot, venting gas. 
© 0000 RAS, MNRAS 000, 000-000 



5 VELOCITY STRUCTURE OF OUTFLOWS 

Figure|7]shows that in all cases the winds have a broad velocity dis- 
tribution extending to > lOOOkms -1 , but most of the wind mass is 
near ~ V c , with relatively little (<C 1% of the mass) at large v > 
500 km s -1 . Observationally, winds in bright ULIRGs have veloci- 
ties similar to those here (Heckman et al. 2000 ; Rupke et al. 2005 ; 
Martin 2005), but in AGN-dominated, late-stage mergers the wind 
velocities typically reach much larger values, ~ 500 — 2000km s _1 
(see c.g.lFeruglio et al.|20T0l|Sturm et al.|2011| [Rupke & Veilleux| 



201l||Greene et aL |20 1 1 1 120 1 2| > | b | This suggests that wind velocity 
may be a useful observational discriminant between starburst and 
AGN driving mechanisms. We find that, without AGN feedback, 
the distribution of wind velocities below the escape velocity but 
above the disk circular velocity is quite flat; above this characteris- 
tic velocity however, it is exponentially decreasing with increasing 

Fig. [8] shows the geometry of this velocity field on differ- 
ent scales, for a specific system (HiZ e) at a given instant (pre- 
coalescence here). On scales comparable to the disks, it is clear that 
there are large non-radial components in the outflows, tracing both 
the orbital motions of the galaxies and the rotation of the disk (with 
local components from e.g. individual star clusters). Since the wind 
is driven from radii throughout the star-forming disk, this rotation 
component is detectable in the wind geometry out to ~ 100 kpc, 
consistent again with what has been seen in ULIRGs ( Martin 2006 ; 
|Gauth ier & Chen 2012). On still larger scales, though, the outflow 
is primarily radial. 

Figs. [9] [10] & [TT] attempt to give an estimate of the ob- 
servable line-of-sight velocity distribution. In Fig. [9] we show the 
mass-weighted line-of-sight velocity distribution of all gas for each 
galaxy, at the snapshot nearest the peak in the starburst (generally 
shortly after nuclear coalescence), after projection onto a random 
axis. We stress that a direct model of e.g. an observed line profile 
requires properly modeling emission and absorption and involves 
three-dimensional line transport, which can differ significantly be- 
tween lines and is outside of the scope of our study here (but see 
e.g. |Cooper et al.|2009| >. But this gives a rough guide to observed 
behaviors, especially if we consider different phases separately. In 
the cold gas especially, we see multiple narrow components with 
separations comparable to the circular velocity; these reflect un- 
relaxed merger kinematics. But in the cold/warm gas, and espe- 
cially in the hot gas, we also see broad wings extending to a few 
hundred km s 



10 



a direct consequence of the winds. In Fig. 
show the same, for just one system but spatially resolved across the 
halo. This makes obvious how the narrow component corresponds 
to different dense filaments and clumpP, but also shows how both 
the narrow offset and broad winds systematically vary across the 
galaxy (owing primarily to projection effects). The broad wings 
have a lower covering factor, but are more obvious out at large 
radii, as material has escaped to larger scales more quickly. Note 
the resemblance between the distributions here and those suggested 
in ob served starbursts (e.g. [ Martin & Bo uche 2009] |Steidel et al.| 
|20101|Gauthier & Chen|2012) . 

In Fig. [TT] we plot the radial outflow velocity of all gas ver- 
sus three-dimensional distance from the galaxy center, at the end 



6 Note that these are the outflows on galactic scales; there are much higher- 
velocity outflows still associated with broad absorption line systems and 
outflows near the AGN itself (for a review, see Veilleux et al. 2005), al- 
though some of these may even have substantial components at ~kpc scales 
{Dunn et al.|2010||Baut"ista et al.|2010) . 
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Figure 7. Top: Thermal bremsstrahlung emission- weighted distribution of gas temperature dLx /dlogTg^, at « 200 Myr after the final coalescence in the 
merger. Different lines in a given panel correspond to gas with different radial velocities (relative to the center of mass). The emission from more rapidly 
outflowing wind material tends to be dominated by the hotter gas with T ~ 10 6-8 K. Middle: Mass- weighted distribution of gas outflow radial velocities 
dm/dv r . The distribution for all gas is divided into the contribution from different phases: cold (T < 2000 K), warm (primarily ionized) (2000 < T < 4 x 10 5 K) 
and hot diffuse (T > 4 x 10 5 K) gas. This outflowing gas consists primarily of a mix of warm and hot gas with some colder material. The warm material 
usually dominates the mass, and is in the form of the filaments and shells seen in Figure[3](the dense "cloudlets" at large radii, while visually prominent, are 
partly numerical and do not dominate by masSJ). Bottom: Mass-weighted density distribution for the wind material, divided into different phases as in the 
middle panel. The cold/warm/hot gas dominates at high/intermediate/low densities, respectively. The "warm" material at very low densities, <C 10 -4 cm -3 is 
previously "hot" material that has adiabatically cooled as it expands; these low densities arise artificially because we do not include an IGM into which the 
wind expands. 



of one example simulation (i.e. at a post-starburst time, although 
the result is qualitatively similar at earlier times). The broad ve- 
locity distribution and mix of orientations/directions is obvious at 
smaller radii, as is the trend towards primarily radial outflow at the 
largest radii (since this is material which has escaped the galaxy). 
A Hubble-like flow develops quickly, and fitting a power-law to 
the median v r versus r gives v r oc r ' 7-1 °. Some of the increase in 
velocity with distance owes to continuous acceleration, but most of 
the trend simply arises because the fastest-moving material escapes 
to the largest radii. We also show the velocity in projection, plotting 
the line-of- sight velocity versus projected distance for all gas. As 
expected, projection effects greatly broaden the distribution. How- 



ever, the trend is similar - fitting a power law |vi os | oc R^ - 10 gives 
a similar scaling. 



6 OUTFLOW RATES AND MASS-LOADING 
EFFICIENCIES 

In Figure [13] we show the wind mass outflow rate during the merg- 
ers. This is typically quantified in terms of its ratio to the star for- 
mation rate M w - m d/M*, i.e. the wind mass per unit mass of stars 
formed. We define the "wind" as material with a positive Bernoulli 
parameter b = (v 2 + 3 c 2 s — v;r sc )/2, i.e., material that would escape 
in the absence of additional forces or cooling. We show the wind 
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Figure 8. Image of the gas velocity field. The gas image from Fig.^is shown, in a face-on (top) and edge-on (bottom) projection, with the (projected in-plane) 
velocity vectors plotted. The vectors interpolate the gas velocities evenly over the image; their length is linearly proportional to the magnitude of the local 
velocity with the longest plotted corresponding to w 500km s _1 . The spatial scale of each image is ±50kpc (left), ±100kpc (middle), ±200 kpc (right). On 
scales near the galaxies the field is complex and the outflow "launching" region traces the entire disk surface, with significant non-radial components from the 
disk and merger orbital motion. On much larger scales the outflow is primarily radial. 



mass above different absolute radial velocity cuts to highlight the 
characteristic velocities 

To first order, the wind mass-loading efficiency does not 
strongly depend on the merger stage. The absolute outflow rate in 
the starburst is very large, ~ 1 — 10 times the SFR for the mod- 
els here, or ^ 10 — 500MQyr -1 in absolute units. But this is not 
proportionally much larger than what is seen for the isolated ver- 
sions of these galaxy models; M W i n d/M* is relatively flat in time. 
In some cases it even drops as the merger begins; this is because 
of the increase in surface densities making the escape of photons 
and hot gas needed to drive the winds less efficient. In Paper III we 
parameterize the dependence of total wind mass-loading on galaxy 
properties (for isolated disks) as 

/Mwind\ „ / Vc(R) \ -(1*0-3) / Sgasfo) \ -0>.5±0.15) 

\ M* / UOOkms-V VlOM pc- 2 / 

(1) 

with a scatter of ~ 50% in normalization. This appears to provide 
a reasonable fit to the values in Fig. [13] 

In the MW case, however, there does appear to be a sharp 
increase in the wind mass-loading at the first-passage and coales- 
cence bursts. This is because, for that model, the gas fraction and so 

7 Note that the first ~ 10 8 yr reflect both out-of-equilibrium initial condi- 
tions and contributions from the "pre-existing" stars in the initial conditions 
to the outflows (not the self-consistently formed stars in the simulation), so 
should be ignored). 
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Figure 9. Line-of-sight velocity distribution of the galaxies, on linear scale. 
Specifically, we project each galaxy along a random axis at a time near the 
peak starburst (just after final coalescence), and take the integrated dra/dvi os 
for all gas, and for the gas separated by temperature (as Fig. [7}. Disturbed 
merger kinematics manifest as asymmetry in the "core," the winds are evi- 
dent in the broad "wings" extending to several hundred km s _1 . We caution 
that this is not the same as an actual observed line profile (there is no ac- 
counting for emission/absorption here), but gives an idea of the line-of-sight 
kinematics in the relevant gas phases. 
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Figure 10. Spatially resolved line-of-sight velocity distributions, for one example simulation (the SMC f merger), at a time just after final coalescence. The 
background image is the simulation gas as Fig.[^ but with inverted brightness (darker is more dense) for clarity. Each box overplots the projected LOSVD of all 
gas as in Fig. [9] averaged over a FWHM 5 kpc Gaussian aperture centered on the box center. The x-axis in each LOSVD box ranges from ±750kms -1 . Dense 
filaments appear as complex multiple narrow features (along with clumps, though we caution these may be numerical artifacts^); lines of sight through larger 
shells have a bimodal appearance; velocity centroid offsets are present at different radii but relatively small (~ 100km s _1 ). The broad velocity distribution is 
evident in the ubiquitous asymmetric tails. 



specific SFR is much lower than any other model; the result shown 
in Paper III is that the "wind" (in the isolated MW case) is mostly 
directly venting hot gas from SNe and stellar winds, rather than any 
entrained material "blown out." It is, in short, below the threshold 
of SFR or "feedback strength" needed to blow out more material. 
But when the specific SFR is boosted in the merger passages, this 
allows it to drive efficient winds. This is not captured in the sim- 



ple scaling for wind mass-loading efficiencies proposed for isolated 
galaxies in Paper III. 

7 EFFECTS OF OUTFLOWS ON STAR FORMATION IN 
MERGERS 

The star formation properties of the simulations here are discussed 
in detail in Paper IV. However, here we wish to consider how these 



© 0000 RAS, MNRAS 000, 000-000 



Starburst Winds & Galaxy Mergers 13 




10 100 1000 

Rproj [kpc] 

Figure 11. Top: Distribution of radial outflow velocities (v r ) for one ex- 
ample simulation (the HiZ e merger, chosen because it features the largest 
velocity range), as a function of three-dimensional distance (r) from the 
galaxy center, at the final time of the simulation. Contours are iso-density 
contours at (0.001, 0.01, 0.1, 0.30.5) times maximum (in the plotted log- 
log coordinates). At smaller radii there is a broad range of outflow speed, 
along with inflow and rotationally supported material, but the median re- 
lation develops a Hubble-like flow. Fitting a power-law, we obtain (v r ) oc 
r o.7-i.o bottom: Same, but in projection. We plot the distribution of the 
absolute value of the line-of-sight velocity versus projected distance from 
the galaxy center (for a random viewing angle). The scatter is much larger 
owing to projection effects, but the scaling is similar |vi os | oc ^pj-L -10 - 



are affected by the outflows generated in the merger. In particular, 
in previous work on galaxy mergers (as discussed in § [TJ it was 
not possible to explicitly resolve wind generation, and instead sim- 
ulations used a variety of sub-grid approaches designed to model 
some "effective" wind scalings. We wish to examine how the con- 
sequences of these winds differ from those here. 

Fig. [T4] shows the SFR versus time in a couple of our mergers, 
compared to simulations with identical initial conditions run using 
a simplified "sub-grid" treatment of stellar winds and the ISM from 
Springel & Hernquist (2003a). In that model, rather than resolve 
the micro- structure of the ISM below ^kpc scales, gas is assigned 
an "effective equation of state" (i.e. effective pressure above a low 
density threshold = 1 cm -3 where the medium is assumed to be- 
come multi-phase) motivated by the interplay of GMC formation 



and destruction, and turbulent driving and heating via stellar feed- 
back. In the "sub-grid" treatment, the lack of resolution of star for- 
mation and feedback means that star formation is assigned statisti- 
cally to gas at much lower densities, with an efficiency that must be 
tuned so that the model lies on the Kennicutt relation. Stellar winds 
are not resolved, so gas particles are instead stochastically kicked 
out of the galaxy at a rate proportional to the SFR, with a fixed ve- 
locity and mass-loading; they are then "free- streamed" (temporar- 
ily decoupled from the hydrodynamics) until they escape to a few 
kpc from the galaxy center (this simply guarantees that the assigned 
mass-loading is the actual wind mass, rather than most of it go- 
ing into e.g. local turbulence). Here, we match the mass-loading to 
the mean measured in each corresponding simulation with the full 
treatment of feedback^] the velocity loading corresponds to a frac- 
tion = 0.1 of SNe energy coupled, but is dynamically irrelevant 
because of the free- streaming condition. For simplicity, we com- 
pare just the Sbc e and MW e models, though the differences are 
robust across other simulations (the differences for the f orbits are 
similar but less pronounced since the orbit produces a less extreme 
starburst). 

We clearly see that the sub-grid wind models "wipe out" some 
of the starbursts at the time of coalescence. This occurs because the 
wind prescription blows material away at all densities, as the gas 
falls into the center; it thus effectively suppresses material actually 
getting into a dense, kpc- scale nucleus in the first place (rather than 
blowing out the material after the starburst begins). It also occurs 
because material cannot be "saved" for the starburst (this requires 
resolved phase structure), nor can dense portions of inflows "resist" 
the outflow occurring simultaneously. A consequence of this is that 
the mergers with certain simplified sub-grid wind models - even 
of gas-rich disks on favorable orbits - have difficulty reproducing 
ULIRGs and other extremely bright merging systems observed (see 
e.g. |Cox et al.|2006b||Dave et al.|2010) . The effect appears more 
severe when gas fractions are large (the low-mass Sbc case), where 
the merger with sub-grid winds here barely produces an enhance- 
ment of the SFR. 

In contrast, in the simulations with resolved feedback the 
cold, dense, compact gas being strongly torqued (and hence tightly 
bound) is difficult to entrain, and so more survives to contribute 
to the starburst, which then efficiently blows out less bound ma- 
terial recycled from inside clouds and/or being contributed from 
more diffuse flows at slightly later times. As shown in Paper IV 
and herein this simultaneously produces large SFR enhancements 
(comparable to "sub-grid" models with zero wind mass-loading) 
and ULIRG/Hyper-LIRG luminosities, while driving winds with 
large mass-loading factors. 

We also see that the post-merger "tail" of star formation is 
significantly different in sub-grid wind models. This is because 
the wind, by design, completely escapes the galaxy; meanwhile 
with no resolved feedback (merely effective pressure) in the disk, 
the material which remains is efficiently exhausted by star forma- 
tion. However, in the simulations with resolved feedback, there is 
a wide range of material in fountains with broad velocity, density, 
and temperature distributions (in addition to the unbound wind); in 
the best-studied systems, this appears to be observationally true as 
well ( Strickland & Heckman 2009 ). Some of this then continuously 
rains down at later times onto the post-merger system, 

Fig.[T5]quantifies this by plotting the distribution of Bernoulli 



8 Though we caution that because of the differing implementations, it is 
difficult to exactly defined "matched" mass-loadings. 
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Figure 12. Galactic super- wind mass outflow rates M w i n d- We compare each disk model (labeled), and both orbits. The mass outflow rate is averaged over 
~2x 10 7 yr intervals, about a dynamical time; wind material is defined by a Bernoulli parameter b > 0. The absolute outflow rates are highest in the HiZ case 
(up to ~ 600 Mq yr -1 over these timescales), but even in the SMC case reach ~ IOMq yr -1 . The peaks broadly follow the starbursts at first couple passages 
and final coalescence. There is surprisingly little orbital dependence in the typical outflow rates. 



parameters b for the gas, at the time just after peak starburst. 
There is a long tail towards b ^> 1 which is the un-bound material, 
but even at large radial velocities v r > lOOkms -1 and distances 
> lOkpc from the disk, a non-trivial fraction of the wind mate- 
rial (tens of percent) is bound (b < 0). This will rain down on the 
disk on a wide range of timescales corresponding to the broad ve- 
locity distribution, from ~ lOOMyr to several Gyr, and re-supplies 
the "tail" of star formation at late times. In contrast to simply con- 
centrating the gas in one place at one time, this makes a complete 
"shutoff" of star formation much more difficult. 

8 DISCUSSION 

In a series of papers, we have implemented detailed, explicit mod- 
els for stellar feedback that can follow the formation and destruc- 
tion of individual GMCs and star clusters: the models include star 
formation only at extremely high densities inside GMCs, and the 
mass, momentum, and energy flux from SNe (Types I & II), stellar 
winds (both "fast" O-star winds and "slow" AGB winds), and radia- 
tion pressure (from UV/optical and multiply- scattered IR photons), 
as well as HII photoionization heating and molecular cooling. Here, 
we extend our models of isolated disks from previous papers to in- 
clude major galaxy mergers. As a first study, we focus on simple, 
global properties, and compare them to those obtained from previ- 
ous generations of simulations which did not follow these processes 
explicitly but instead adopted a simplified "effective equation of 
state" sub-grid model of the ISM. 

With explicit feedback models, super-winds are generated in 
all passages with outflow rates from ^ 10 — 500 Mq yr -1 . The sim- 
ulated outflow rates are suggestively similar to observations of star- 



burst winds in many observed ULIRGs and merging galaxies (com- 



pare e.g.|Heckman et al.|1990 2000l|Rupke et al.|2005||Coil et al.| 


2011 


Grimes et al. 2009, Rubin et al. 2011, Martin 2005 2006 


1999 


|Sato et al.|2009||Gauthier & Chen|2012|). 



Although the absolute outflow rates can be enormous in the 
starburst, we find that the mass-loading efficiency - i.e. the outflow 
rate per unit star formation rate - is broadly similar for each galaxy 
model to the isolated version of that system. In other words, out- 
flow efficiencies - on average - appear to depend more strongly on 
global dynamical properties (escape velocity, effective radius) than 
on the instantaneous dynamical state of the system. The character- 
istic mass-loading efficiency scales inversely with the escape veloc- 
ity of the systems, increasing from about unity in massive systems 
(~ 10 11 Mq) to ~ 10 — 20 in SMC-mass dwarfs. However since 
the absolute SFRs are much larger in massive systems, the absolute 
outflow rates tend to increase with galaxy mass (for gas-rich sys- 
tems). We caution that, as for isolated systems in Paper III, there is 
very large variability in the instantaneous mass-loading efficiency 
M W i n d/M* - at least an order of magnitude scatter. 

We make predictions for the distribution of velocities, densi- 
ties, and temperatures/phases of the outflows, which extend those 
shown in Paper III for isolated disks. The distribution of veloci- 
ties is broad, as shown for isolated systems in Paper III, but ex- 
tends characteristically to the escape velocity (a couple times the 
maximum circular velocity), with a long (but relatively low-mass) 
tail of material in higher- velocity components. This is also similar 
to observations - typical velocities are a couple to a few hundred 
kms -1 (references above). The winds are characteristically multi- 
phase, with a large fraction of the mass in each of the cold, warm 
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Figure 13. Average galactic wind mass-loading efficiency (= M^- m ^/M nQ ^; where M w j nc i = f M w i n d and M new = / M*) for each galaxy model (f/e orbits 
shown as solid/dotted lines as Fig.[T2). As shown in Paper III, the mass-loading increases going from high-mass (HiZ, MW) systems to lower-mass (Sbc) and 
dwarf (SMC) galaxies. The mechanisms dominating the outflow also transition from radiation pressure to SNe, respectively. As isolated disks, the HiZ, Sbc, 
MW, SMC models have M w ind ~ (1,3 — 5,1 — 3, 10 — 20) M* . Similar values appear here, with weak time-dependence. Although the absolute outflow rates 
in starbursts in Fig.[^|are very large, they follow from large SFRs; the mass loading efficiency is more sensitive to global galaxy properties than merger stage. 
We also consider the mass-loading of material with various cuts in the outflow radial velocity v r (in kms -1 ); the distribution of outflow velocities is very 
broadP 



ionized, and hot (pressure-supported) phase, and a broad range of 
densities spanning several orders of magnitude. The cold phases are 
predominantly accelerated by radiation pressure, and become more 
prominent in the more gas-rich, higher-density systems, including 
the central regions of the merger, while gas heated to high tempera- 
tures by SNe ejecta is more prominent in low-mass and/or gas-poor 
systems, and in the more diffuse (volume-filling) outflow. 

In contrast to isolated disks, the wind kinematics are (un- 
surprisingly) complex. The winds do not entirely originate in the 
nuclear starburst, but over the entire surface of the disks as the 
merger proceeds; not only do the extended disks contain non-trivial 
star formation (see Paper IV), but the material there (being lower- 
density, less tightly-bound, and less strongly-torqued) is easier to 
accelerate into the wind (or entrain in the outflow emerging from 
the very center). As such, near the disks the kinematics retain mem- 
ory of the orbit and disk rotation (as well as more chaotic compo- 
nents from individual super- star clusters on sub-kpc scales), as has 
been mapped in detailed studies of local merging systems ( |Mar-| 
tin 2006 ). On larger scales the outflow becomes primarily radial, 
but tends towards a unipolar or bipolar structure; however since the 
disk orientations are changing during the merger, this leads to mul- 
tiple overlapping shells/bubbles with different orientations at dif- 
ferent radii, also similar to observations (Gauthier & Chen 2012| ). 
Because the winds contain material at a wide range of velocities, 
different outflows "catch up" to one another, and a complex three- 
dimensional structure develops. This also imprints a Hubble flow- 
like pattern of larger velocity material at larger radii (even without 



the effects, present here, of continuous thermal and radiation pres- 
sure acceleration), similar to that observed (references above and 
|Martin & Bouche|2009||Steidel et al.|2010| ). All gas phases include 
components extending to velocities of several hundred kms -1 , but 
the cold/warm phases also show narrow features with smaller off- 
sets ~ 100 kms -1 that indicate both merger kinematics and the 
structure of individual streams and entrained/accelerated shells of 
cool gas (though the breakup of these shells into some small cold 
clumps at large radii from the galaxy is subject to significant numer- 
ical caveats^), while the hot components show a smoother, broader 
velocity distribution. 

This multi-phase nature of the wind and its driving across the 
disk are critical to the fact that systems can simultaneously drive 
winds with large average mass-loading factors and also avoid "wip- 
ing out" all structure inside the disk - including the kpc-scale gas 
concentrations that power the starburst itself! In contrast, sub-grid 
models which do not resolve the generation of winds but simply in- 
sert some mass-loading by hand can produce very different results. 
We show that some implementations of these subgrid models can 
suppress the merger-induced enhancement of the SFR - which can 
(when mass-loadings are large) make it difficult to form ULIRGs or 
bright sub-millimeter galaxies. This may be related to some known 
difficulties reproducing these rare populations in cosmological sim- 
ulations (e.g. |Dave et al.|2010| |Hayward et al.|20lTa] >. It is clear 
that careful treatment of subgrid wind models is necessary if one 
wishes to properly resolve the dynamics of gas flows and SFR en- 
hancements within galaxies. 
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Figure 14. Star formation rate versus time for our simulations Sbc and MW 
e ("full feedback"; for the remaining simulations, see Paper IV Fig. 8). We 
compare to simulations run with identical initial conditions but a sub-grid 
model for both the ISM phase structure ("effective equation of state") and 
winds (particles "kicked out" of the galaxy at fixed mass-loading, matched 
to the mean value in the corresponding "full feedback" simulation). The 
sub-grid wind model "wipes out" substructure, blowing out gas before it 
collects in the central kpc and dramatically suppresses the multi-phase star- 
bursts. Also, because gas in the winds is completely ejected from the galaxy, 
the post-merger SFR decays too quickly in the sub-grid model. The differ- 
ence is larger in the Sbc case because the sub-grid wind is assigned a larger 
mass-loading (= 5, vs. = 1 for the MW model), but evident in both. 



It is also important that a significant fraction of the wind ma- 
terial is not unbound, and falls back into the disk over the cou- 
ple Gyr after the final coalescence starburst, leading to a slow, ex- 
tended "tail" in the starburst decay. Especially in prograde mergers, 
this can greatly enhance the magnitude and duration of post-merger 
star formation, relative to the older models which treat the ISM and 
feedback physics in a "sub-grid" manner. This has important im- 
plications for "quenching" of star formation in massive galaxies. If 
quenching were possible without the presence of some additional 
feedback source - say, from an AGN - then the simulations here 
are the most optimal case for this. They are cosmologically iso- 
lated galaxies, so there is zero new accretion; moreover, an equal- 
mass merger represents the most efficient means to exhaust a large 
amount of gas quickly via star formation, much moreso than an iso- 
lated disk (see e.g. |Hopkins et al.|2008c|a| ). But we find that with 
the presence of stellar feedback, many of our models - including 
the already gas-poor MW-like system, maintain post-merger SFRs 
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Figure 15. Distribution of Bernoulli parameters (binding energy) of the 
gas just after the peak of the merger-induced starburst. For each model (f 
merger shown, but e is very similar) we plot the distribution of mass per unit 



b I Vg SC , where b = ( v 2 + 3 cj 



c )/2, so that b > corresponds to material 



which would escape in the absence of additional pressure forces and cool- 
ing. Black lines compare all gas, outside of different radii from the center 
of mass. Red lines correspond to the same limiting radii but restricted to gas 
with a radial outflow velocity v r > 100km s -1 . There is clearly a large tail 
of b > un-bound gas in each case; the secondary peak at b < is virial- 
ized gas. But even for gas with large outflow speeds > 100km s _1 and at 
radii > 10 kpc, there is a non-negligible (~ 10 — 30%) fraction with b < 0, 
which will fall back into the disk. 



nearly as large as their steady-state pre-merger SFR. The systems 
simulated here would take several Gyr to cross the "green valley" 
and turn red, much longer than the <Gyr quenching timescale re- 
quired by observations (see M artin et aT| 2007 ; Snyder e t al.|20lf] ). 
Far from resolving this by gas expulsion, stellar feedback makes the 
"quenching problem" harder. As shown in Moster et al. (201 1), ad- 
dition of realistic gas halos around the merging galaxies (even with- 
out continuous accretion) only further enhances the post-merger 
SFR. This is the short-timescale manifestation of a general prob- 
lem in cosmological simulations; over a Hubble time, recycled ma- 
terial from galaxy progenitors is re-captured, and leads to large star 
formation rates and excessive stellar masses in high-mass galaxies 
(Oppenheimer et al. 2010). If gas is to be swept out of galaxies 
efficiently after a merger or starburst, the models imply that some 
other form of feedback - perhaps from bright quasars - is neces- 
sary. This is also suggested by observations of late-stage mergers, 
which find that in the AGN-dominated systems at quasar luminosi- 
ties, outflow masses are enhanced and the outflow velocities reach 
~ 1000km s _1 , larger than those we find driven by stellar feedback 
(|Feruglio et al.|2010[|Tremonti et al.|2007l|Sturm et al.|201 l[|Rupke| 
|& Veilleux|2011| ). These high velocities (well above the escape ve- 
locity) may provide a unique signature of AGN-driven outflows, 
since our simulations drive very little mass to such high values. But 
we caution that the relative timing of AGN and starbursts is uncer- 
tain (though the simulations suggest AGN follow the starburst; see 
Hopkins 2011 ), and AGN (owing to their complicated duty cycles) 
may not still be active when AGN-driven winds are observable. 

In a companion paper, we will examine the star clusters 
formed in these simulations. The mass/luminosity distribution, spa- 
tial locations, formation time distribution, and physical properties 
of these clusters represent a powerful constraint on small-scale 
models of the ISM and star formation. 

We have also restricted our focus to major mergers. Studies 
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of mergers with varying mass ratios suggest that the qualitative be- 
haviors discussed here should not depend on mass ratio for ratios 
to about 3:1 or 4:1, and even at lower mass ratios they can be con- 
sidered similar but with an "efficiency" of inducing starbursts and 
violent relaxation that scales approximately linearly with mass ra- 
tio (Hernquist & Mihos|| 1995] |Cox et"al][2008l [Naab & Burkert] 
2003 ; Younger et al. 2008a). Since the simplest properties of the 
winds (e.g. their mass loading and characteristic velocities) seem 
to scale relatively simply between isolated disks and major merg- 
ers, we expect that minor mergers will represent an intermediate 
case. 

We note that recent studies comparing cosmological simula- 
tions done with GADGET and the new moving mesh code AREPO 
(Springel 2010) have highlighted discrepancies between grid codes 
and smoothed particle hydrodynamics (SPH) in some problems re- 
lated to galaxy formation in a cosmological context (Agertz et al. 
Vogelsberger et aL|[20T2| |Sijacki et aL]|2012[ |Keres et al. 



2007 



2012 



Bauer & Springel 2012] |Torrey et al.||201l| ). However, we 



have also performed idealized simulations of mergers between in- 
dividual galaxies and found excellent agreement between GADGET 
and AREPO for e.g. gas-inflow rates, star formation histories, and 
the mass in the ensuing starbursts (Hayward et al.|20lTb] >. More- 
over, in [Hopkins ( 2012a) we show that many of these discrep- 
ancies can be resolved with small modifications to the equations 
of motion, and test the differences in galaxies with the full feed- 
back models presented here. Although subtle numerical issues can 
influence quantities like fluid mixing and hence hot gas cooling, 
we show that the SFRs agree very well and wind masses agree to 
within a factor of two, much smaller than the differences if we re- 
move feedback. Finally, the differences in numerical methods are 
also minimized when the flows of interest are supersonic (as op- 
posed to sub-sonic), which is very much the case here (Kitsionas 
|et al.|20091|Price & Federrath|20101 [Bauer & Springel|2012| ). 

These new models allow us to follow the structure of the gas 
in the central regions of starburst systems at high resolution. This 
makes them an ideal laboratory to study feedback physics under 
extreme conditions in, say, the center of Arp 220 and other very 
dense galaxies. We have also, for clarity, neglected AGN feedback 
in these models, but we expect it may have a very significant effect 
on the systems after the final coalescence (e.g. |Di Matteo et al.| 
|2005[ |Springel et"aL]|2005a|b) . With high-resolution models that 
include the phase structure of the ISM, it becomes meaningful to 
include much more explicit physical models for AGN feedback. 
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Figure Al. Morphology of the gas in the isolated (non-merging) SMC 
model, taken from Hopkins (2012a). Colors encode gas temperature as 
Figs.[^&[3]in face-on {top) and edge-on {bottom) projections. Left: The sim- 
ulation uses the standard "density-entropy" formulation of the SPH equa- 
tions of motion, our "default" implementation in this paper. Right: The 
simulation uses the alternative "pressure-entropy" formulation developed 
in Hopkins (2012a), which resolves some known problems in the "stan- 
dard" SPH method treating fluid mixing instabilities and contact discon- 
tinuities (e.g. the Kelvin-Helmholtz instability), but is otherwise identical. 
The two are largely similar - the pressure-entropy formulation produces 
slightly less-sharp transitions between molecular features and hot bubbles 
(owing to mixing), but the difference is small. 
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Figure A2. Star formation rate and stellar wind mass-loading (as Fig.[T3), 
for the same isolated SMC models in Fig. |A1| again taken from Hopkins 
{2012a}. The time-averaged SFR differs by only ~ 20% in the two SPH im- 
plementations (time variations can be larger, but these are stochastic). The 
mass-loading is lower in the pressure-entropy case by a factor ~ 1.5 — 2.0, 
because the increased mixing along phase boundaries brings more (initially) 
cold, dense gas into "hot" gas, enhancing its cooling rates (before it can ac- 
celerate into an outflow). The time-dependence and qualitative behavior is 
similar. 
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APPENDIX A: NUMERICAL TESTS OF THE SPH 
METHOD AND WIND PHASE STRUCTURE 

In this appendix, we examine the robustness of the numerical meth- 
ods used here - in particular, we wish to study how the small-scale 
phase structure of the outflows can be affected by detailed of the 
methodology. 

Our default simulations in this paper use the standard 
"density-entropy" formulation of the SPH equations of motion 
in GADGET from |Springel & Hernquist| p052l >. This formula- 
tion manifestly conserves momentum, energy, angular momentum, 
and entropy (in the absence of sources/sinks), and has a number 
of additional advantages, but produces a resolution- scale "surface 
tension"-like error term at contact discontinuities, which has the 
effect of suppressing the growth of some fluid mixing instabilities, 
and has been the subject of much discussion in the literature (see 
Agertz et al. 2007 ; Read & Hay field 2012 and references therein). 

Since the multi-phase winds may well be subject to exactly 
these instabilities, we have re -run a subset of our simulations using 
the newer "pressure-entropy" SPH formulation described in |Hop-| 
|kins| ( [2012a| ), which is shown there to give dramatically improved 
results on in situations with fluid mixing around contact discon- 
tinuities (e.g. the Kelvin-Helmholtz and Rayleigh-Taylor instabili- 
ties) while retaining excellent conservation properties, and includes 
a number of additional improvements to the treatment of artficial 
viscosity (see Cullen & Dehnen 2010), SPH smoothing kernel ac- 
curacy (Dehnen & Aly 2012] ), and timestep communication rele- 



Figure A3. Morphology of the gas in the merging SMC f model, with col- 
ors encoding temperature as Figs. ^[3] & |A1| As Fig. |Al| we compare two 
projections of the same simulation (at an identical time near coalescence). 
The "density-entropy" case is our standard model in the text. The pressure- 
entropy case has been re-run with the revised pressure-entropy SPH for- 
mation, as well as an improved SPH smoothing kernel, artificial viscosity 
scheme, and includes UV background heating and self-shielding. The mor- 
phology of the hot gas shells/bubbles, warm outflowing filaments/shells at 
~ 10 — 100 kpc and tidal tails is very similar. However, the previously cool 
gas at large radii C$> 10 kpc) is now photo-heated to "warm" temperatures 
(in e.g. the tidal tails). The biggest difference is that the small, cold clumps 
at the largest radii (which were never self-gravitating and appeared to form 
in the outflow in the density-entropy models) do not appear in the new code. 
They are a numerical artifact of the "density-entropy" equation of motion. 
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Figure A4. Phase structure and density distribution of the gas in different components in the simulations in Figs. |Al|A3] in the style of Fig. [5] Top: Mass- 
weighted density PDF for all gas in the simulation, and the gas within the multi-phase star-forming disk, the extended ionized disks and halo, and wind/outflow 
material. Middle: Mass- weighted density PDF within the star forming disk, divided into the cold/warm/hot phases. Bottom: Volume- weighted density PDF 
in the star-forming disk. We compare both the isolated SMC progenitor disks from Fig. |A1| and the merging SMC f models in Fig. |A3| each run with our 
"default" density-entropy formulation from the text and re-run with the improved pressure-entropy formulation from Hopkins ( 2012a). As suggested from 
the Figures, the pressure-entropy formulation leads to more mixing with reduces the "hot" wind mass by a factor ~ 2. The phase distribution within the star 
forming disk is nearly identical (down to < 1 % within the tails of the distributions) - we stress that the elimination of the cold "blobs" at large radii in the 
outflow in Fig. |A3| does not apply to the GMCs and star clusters forming within the disk (which are well-resolved and self-gravitating, unlike the cold blobs 
in outflow). Note that the merger run with the pressure-entropy formulation self-consistently includes UV background heating, so the cold gas distribution 
self-consistently "cuts off" at a density < 0.01 cm -3 , very similar to where our post-processing estimates truncate the distribution in the text. Differences in 
the artificial viscosity and other numerical parameters appear to have small effects here. 



vant for treating extremely high Mach-number shocks ( Sai toh &] 
|Makino|2009l|Durier & Dalla Vecchia|2012| >. For extensive numer- 
ical tests demonstrating accurate treatment of these instabilities, see 
|Hopkins| ( [20T2aT ). 

To test whether these subtleties may be strongly influencing 
our results, we first consider an isolated, star-forming disk (the 
progenitors in the SMC model mergers), which was analyzed in 
Hopkins ( 2012a). The simulations are otherwise identical (and in- 
clude the identical physical prescriptions to the merger simula- 
tions herein), but adopt either the density-entropy "standard" SPH 
or newer Hopkins (2012a) pressure-entropy form (here we keep 
the kernel and all other properties fixed between the simulations). 
Fig. |A1 1 shows that the morphologies of the isolated disks are very 
similar; there is a small difference where the pressure-entropy for- 
mulation leads to increased mixing along phase boundaries owing 
to the instabilities above, producing less-sharp divisions between 
molecular regions and hot bubbles, and the star-forming and "ex- 
tended" disks. Fig. |A2| shows that the SFR differs very little - 
the time-averaged SFR differs by only ~ 20% (larger variations 
at any given time appear, but these are consistent with stochastic 
variation that occur re-running simulations with the same numeri- 
cal method). The total wind mass-loading is, however, significantly 
affected, and is lower in the pressure-entropy run by a factor of 



~ 1.5 — 2. This is because the increased mixing adds some cold 
gas to the hot medium (rather than keeping it separate in small, 
cold clumps), which then substantially increases the cooling rate of 
the hot gas. So we caution that the absolute wind mass-loading fac- 
tors should be considered uncertain at the factor ~ 2 level; but this 
is comparable to the level of uncertainty that stems from different 
numerical implementations of the same feedback mechanisms, as 
shown in Hopkins et al. ( 2 012d| >. 

We next compare a merger of these galaxies, our SMC f 
model. Fig. | A3 1 shows the morphology of the winds at large radii in 
the merger (as Fig. [3}, for our standard "density-entropy" SPH for- 
mulation, and the newer "pressure-entropy" formulation. In the lat- 
ter case, we now use the fully updated code from Hopkins ( 2 012a| ), 
with a more accurate artificial viscosity scheme, SPH kernel, and 
timestep limiter. In the newest version of the code, we have also 
implemented heating from the z — UV background as tabulated 
in Faucher-Giguere et al. ( 2008), accounting for self- shielding, and 
we include this as well since it may well alter the phase distribution 
of the gas at large radii. The resulting differences at large radii are 
much more visually striking than those in the star-forming disk: the 
cold clumps or "blobs" at large radii disappear in the new simula- 
tion. We stress that the morphology of the smooth gas - both the 
diffuse volume-filling hot gas and also the warm shells/filaments 
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Figure A5. Phase structure, velocity, and density distributions of the winds of the same simulations as Fig. |A4| in the style of Fig. [7] Top: Mass-weighted 
distribution of gas outflow radial velocities, for all gas and divided into cold/warm/hot phases. Bottom: Mass-weighted density distribution of material in the 
wind, divided by phases. Again we compare the isolated SMC model and SMC f merger, at the same instant in time, but with either our default "density- 
entropy" SPH or the revised "pressure-entropy" SPH. The resulting differences in the winds are larger than in the star-forming disks in Fig. |A4| However 
they are still qualitatively similar. The temperature and density distributions have similar peaks and widths, but differ in their tails. As expected from the 
morphologies in Figs. |Al| & |A3| the cold "blobs" at large radii are more efficiently mixed with the warm outflow in the pressure-entropy formulation, leading 
to a smaller cold gas contribution in the wind (but this is not dominant, in any case). In the merger simulation, this mixing also enhances the cooling of the hot 
gas, so it is reduced as well (making the warm gas component even more dominant). In the merger, the total outflow mass at this stage is systematically lower 
in the pressure-entropy formulation by about ~ 40%. 



- is nearly identical. And within the star-forming disk, GMCs 
still form in very similar fashion. It is largely these cold blobs at 
large radii that are altered. These are not self-gravitating, and are 
marginally resolved. They form by "breakup" of filamentary struc- 
tures in the wind, in the density-entropy formulation. But this ap- 
pears to be a numerical artifact of the density-entropy formulation, 
specifically the "surface tension"-like term (which drives a bead- 
ing effect along accelerating filaments that leads to breakup). The 
ionizing background suppresses their cooling to cold temperatures 
(though it has little effect on the morphology). 

These distinctions are more evident in the distribution of gas 
phases. In Fig. |A4[ we repeat Fig. [5] and examine the density dis- 
tribution of gas in the winds, ionized disk, and star-forming disk, 
as well as the specific phase breakdown within the star-forming 
disk. We compare both the isolated and merging SMC models with 
both density-entropy and pressure-entropy SPH. First, we empha- 
size that all of our qualitative conclusions appear robust. There are 
some quantitative changes, but most of these are in the tails of the 
density distributions, relevant at the sub-percent level. Within the 
star-forming disk, we see that the phase distributions are very sim- 
ilar in both implementations of SPH. Furthermore, including the 
ionizing background self-consistently has very little effect, except 
to truncate the cold gas density distribution at just about the den- 
sity where our previous simple post-processing estimate led us to 
truncate the distributions. The predicted properties within the disk 
appear very robust. 



In Fig. |A5[ we extend this comparison to the velocity and den- 
sity distributions of the wind material (as in Fig. [TJ. Here, we see 
larger differences, as expected from the previous Figures. In the 
pressure-entropy SPH, enhanced mixing of cold and hot phases de- 
creases the relative importance of both in the wind, and enhances 
the relative importance of the warm-phase gas. The presence of an 
ionizing background also contributes to this. However, in each case 
we stress that warm-phase gas already dominated the outflow, so 
this conclusion is robust. The velocity distributions are altered, but 
only at a modest level - there is still a broad velocity distribution in 
all phases (in fact, in the models here at the specific time analyzed, 
there may be somewhat more cold material at very large velocities, 
even though there is less overall in the outflow). 

We have specifically chosen to focus on the SMC case here, 
because its outflow, being predominantly "hot phase" gas (but fea- 
turing some cold blobs in the density-entropy runs) is most likely 
to be strongly affected by the detailed of the numerical method 
and fluid mixing. We have, however, also re-run lower-resolution 
versions of the HiZ e and MW e mergers with both the density- 
entropy and pressure-entropy formulations of SPH (and have ex- 
perimented with a wide range of artificial viscosity and smoothing 
kernel implementations, as well as ionizing background strengths, 
in the isolated disk progenitors of each). In all these cases, we find 
the sense of the difference between SPH formulations is identical 
to that described above - however, the quantitative magnitude of 
the difference is smaller in each case (see Fig. |A6| ). The total wind 
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Figure A6. Wind thermal+metal line emission morphology, as Fig. [4] but 
for the SMC f merger re-run with the revised "pressure-entropy" formu- 
lation of SPH (bottom), as well as a HiZ e merger run with the pressure- 
entropy SPH but lower resolution (lOx lower particle number). The large- 
scale wind behavior - at least in these simulations which do not include 
a cosmological IGM - is not strongly sensitive to the numerical method 
(modulo the small clumps in Fig. [A3}. 

mass-loading, for example, is only lower by ~ 20% in the MW run, 
and actually appears to be slightly higher in the HiZ case. 

To summarize, the qualitative conclusions and results pre- 
sented in the main text all appear robust to the details of the numer- 
ical method, except for the presence of cold "blobs" which form in 
the outflow. These features are likely driven by the same numerical 
artifacts that have already been identified as causing the breakup 
of inflows in cosmological simulations into similar "clumpy" mor- 
phologies (see Keres et al. 2012; Nelson et al. 2013 ). We stress, 
however, that these do not dominate the outflows by mass. More- 
over, they do not appear until the outflow has essentially escaped 
the disk. At this point, as we have emphasized in the text, the de- 
tailed phase structure predicted should not be taken too seriously, 
since there is no IGM into which the winds expand. However, this 
should strongly caution against over-interpretation of the detailed 
phase distribution in cosmological simulations with outflows (many 
of which show similar features). 

Quantitatively, differences in the numerical method and ioniz- 
ing background can lead to systematic changes in our predictions 
for basic wind properties (mass-loading factors, total mass in dif- 
ferent phases, maximum velocities) at the factor of ^ 2 level. This 
is significant. As shown in Paper III, similar uncertainties can arise 
owing to the manner in which stellar feedback physics is imple- 
mented. And this is also comparable to genuine physical uncertain- 
ties, for example, the strength of feedback can vary at this level 
owing to plausible variation in the stellar initial mass function, or 
the presence of an equipartition magnetic field and/or cosmic rays 



(see e.g. |Uhlig et al.|2012l|Pakmor & Spri ngel 2012 and references 
therein). So we strongly emphasize that considerably more sophis- 
ticated simulations are needed (along with improved observational 
constraints) before any "precision" predictions with accuracy much 
better than a factor of ^ 2 can be made. 
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